## Work-from-home, electricity, and water: Evidence from COVID-19 in Qatar
## Author: David Bernstein 

## Clear global environment 
rm(list=ls())

## Load in packages - using pacman for easy loading 
pacman::p_load(readxl,matrixStats,reshape2,Hmisc,dplyr,stargazer,plm,tseries,quantreg)

## define dummy variable (COVID or STR) and use of quantile or fe regressions
USE_DUMMY <- TRUE 
QUANTILE  <- TRUE

## load data 
load("data.RData")

## Basic panel regression analysis
## deal with potential zeros -- make zeros 1 in order to pass to log
elec_comm_panel[which(elec_comm_panel$CONSUMPTION==0),       "CONSUMPTION"]       <- 1
elec_resi_panel[which(elec_resi_panel$CONSUMPTION==0),       "CONSUMPTION"]       <- 1
water_comm_panel[which(water_comm_panel$CONSUMPTION==0),     "CONSUMPTION"]       <- 1
water_resi_panel[which(water_resi_panel$CONSUMPTION==0),     "CONSUMPTION"]       <- 1

## Use for log transformation of STR 
# elec_comm_panel[which(elec_comm_panel$stringency_index==0),  "stringency_index"]  <- 1
# elec_resi_panel[which(elec_resi_panel$stringency_index==0),  "stringency_index"]  <- 1
# water_comm_panel[which(water_comm_panel$stringency_index==0),"stringency_index"]  <- 1
# water_resi_panel[which(water_resi_panel$stringency_index==0),"stringency_index"]  <- 1

elec_comm_panel[, "t"]  <- as.numeric(elec_comm_panel[,  "DATE"])
elec_resi_panel[, "t"]  <- as.numeric(elec_resi_panel[,  "DATE"])
water_comm_panel[,"t"]  <- as.numeric(water_comm_panel[, "DATE"])
water_resi_panel[,"t"]  <- as.numeric(water_resi_panel[, "DATE"])

if (USE_DUMMY == TRUE ){
      EQUATION_C <- log(CONSUMPTION) ~ COVID_DUMMY             + factor(MONTH)
      ROW_LABS <- c("COVID","Aug","Jul", "Jun","Mar","May","Oct","Sep","t","constant")
}else{EQUATION_C <- log(CONSUMPTION) ~ I(stringency_index/100) + factor(MONTH) 
ROW_LABS <- c("STR","Aug","Jul", "Jun","Mar","May","Oct","Sep","t","constant")}

## Causes many convergence issues
# if (USE_DUMMY == TRUE ){
#   EQUATION_C <- log(CONSUMPTION) ~       COVID_DUMMY           + factor(MONTH) + factor(YEAR)
# }else{EQUATION_C <- log(CONSUMPTION) ~ I(stringency_index/100) + factor(MONTH) + factor(YEAR)}

## Panel estimators
plm_fe_e_c  <- plm(EQUATION_C , elec_comm_panel, effect = "individual", model = "within")
plm_fe_w_c  <- plm(EQUATION_C,  water_comm_panel,effect = "individual", model = "within")

if(QUANTILE == TRUE){
  ## Quantile Estimator
  rq_fe_e_c  <- rq(EQUATION_C , data = elec_comm_panel,  tau = 1:9/10)
  rq_fe_w_c  <- rq(EQUATION_C,  data = water_comm_panel, tau = 1:9/10)
  
  pdf(file="rq_elec_c_str.pdf")
  plot(summary(rq_fe_e_c, se = "iid"))
  dev.off()
  
  pdf(file="rq_water_c_str.pdf")
  plot(summary(rq_fe_w_c, se = "iid"))
  dev.off()}else{print("NO QUANTILE")}

## Panel Regressions for subsets of the data for Water and Electric
## Take subsets of Electricity 
## Qatari   
elec_resi_panel_qat_all  <- subset(elec_resi_panel, 
                                   (Customer.Type=="Qatari Owners") & 
                                     (Premises.Category== "Villa (Residential)" |
                                        Premises.Category== "Flat (Residential)")  )

elec_resi_panel_qat_vil  <- subset(elec_resi_panel, 
                                   (Customer.Type=="Qatari Owners") & 
                                     Premises.Category== "Villa (Residential)"  )

elec_resi_panel_qat_flat <- subset(elec_resi_panel, 
                                   (Customer.Type=="Qatari Owners") & 
                                     Premises.Category== "Flat (Residential)"   )  

## Regular Customers 
elec_resi_panel_reg_all  <- subset(elec_resi_panel, 
                                   (Customer.Type=="Regular Customers") & 
                                     (Premises.Category== "Villa (Residential)" |
                                        Premises.Category== "Flat (Residential)")  )

elec_resi_panel_reg_vil  <- subset(elec_resi_panel, 
                                   (Customer.Type=="Regular Customers") & 
                                     Premises.Category== "Villa (Residential)"  )

elec_resi_panel_reg_flat <- subset(elec_resi_panel, 
                                   (Customer.Type=="Regular Customers") & 
                                     Premises.Category== "Flat (Residential)"   )
## Take subsets of Water 
## Qatari  
water_resi_panel_qat_all  <- subset(water_resi_panel, 
                                    (Customer.Type=="Qatari Owners") & 
                                      (Premises.Category== "Villa (Residential)" |
                                         Premises.Category== "Flat (Residential)")  )

water_resi_panel_qat_vil  <- subset(water_resi_panel, 
                                    (Customer.Type=="Qatari Owners") & 
                                      Premises.Category== "Villa (Residential)"  )

water_resi_panel_qat_flat <- subset(water_resi_panel, 
                                    (Customer.Type=="Qatari Owners") & 
                                      Premises.Category== "Flat (Residential)"   )  

## Regular Customers 
water_resi_panel_reg_all  <- subset(water_resi_panel, 
                                    (Customer.Type=="Regular Customers") & 
                                      (Premises.Category== "Villa (Residential)" |
                                         Premises.Category== "Flat (Residential)")  )

water_resi_panel_reg_vil  <- subset(water_resi_panel, 
                                    (Customer.Type=="Regular Customers") & 
                                      Premises.Category== "Villa (Residential)"  )

water_resi_panel_reg_flat <- subset(water_resi_panel, 
                                    (Customer.Type=="Regular Customers") & 
                                      Premises.Category== "Flat (Residential)"   )

subset_names <- c("elec_resi_panel_qat_vil", "elec_resi_panel_qat_flat",
                  "elec_resi_panel_reg_vil", "elec_resi_panel_reg_flat",
                  "water_resi_panel_qat_vil", "water_resi_panel_qat_flat",
                  "water_resi_panel_reg_vil", "water_resi_panel_reg_flat")

## Tau = 0.50 and Mean from OLS  
for (i in 1:8){
  print(paste(i," ",subset_names[i]," ",sep=""))
  DATA   <- get(subset_names[i])
  STORE  <- plm(EQUATION_C ,data = DATA , effect = "individual", model = "within")
  STORE0 <- lm(EQUATION_C ,data = DATA)    
  assign(paste("plm_fe_sub_",i,sep = ""), STORE)
  assign(paste("lm_fe_",i,sep = ""),     STORE0)
  
  if(QUANTILE == TRUE){
    STORE1 <- rq( EQUATION_C ,data = DATA , tau = 1:9/10)   
    assign(paste("rq_fe_sub_", i,sep = ""), STORE1)
    
    pdf(file= paste("rq_fe_sub_str_", i,".pdf" ,sep = ""))
    plot(summary(STORE1, se = "iid"))
    dev.off()
    rm(STORE1)}else{print("NO QUANTILE")}
  rm(STORE, DATA)  
}

plm_fe_e_c  <- plm(EQUATION_C, elec_comm_panel, effect = "individual", model = "within")
plm_fe_w_c  <- plm(EQUATION_C, water_comm_panel,effect = "individual", model = "within")
lm_fe_e_c   <- lm(EQUATION_C, elec_comm_panel)
lm_fe_w_c   <- lm(EQUATION_C, water_comm_panel)

stargazer(plm_fe_sub_1,plm_fe_sub_2,plm_fe_sub_3,plm_fe_sub_4, plm_fe_e_c, covariate.labels = ROW_LABS)
stargazer(plm_fe_sub_5,plm_fe_sub_6,plm_fe_sub_7,plm_fe_sub_8, plm_fe_w_c, covariate.labels = ROW_LABS)

stargazer(lm_fe_1,lm_fe_2,lm_fe_3,lm_fe_4, lm_fe_e_c, covariate.labels = ROW_LABS)
stargazer(lm_fe_5,lm_fe_6,lm_fe_7,lm_fe_8, lm_fe_w_c, covariate.labels = ROW_LABS)


## Tau = 0.50
tau <- 0.50

for (i in 1:8){
  print(paste(i," ",subset_names[i]," ",sep=""))
  DATA   <- get(subset_names[i])
  
  if(QUANTILE == TRUE){
    STORE1 <- rq( EQUATION_C ,data = DATA , tau = tau)
    assign(paste("rq_tau_sub_", i,sep = ""), STORE1)
    
    rm(STORE1)}else{print("NO QUANTILE")}
  
  rm(DATA)  
}

rq_tau_e_c  <- rq(EQUATION_C , data = elec_comm_panel,  tau = tau)
rq_tau_w_c  <- rq(EQUATION_C,  data = water_comm_panel, tau = tau)

stargazer(rq_tau_sub_1,rq_tau_sub_2,rq_tau_sub_3,rq_tau_sub_4,  rq_tau_e_c, covariate.labels = ROW_LABS)
stargazer(rq_tau_sub_5,rq_tau_sub_6,rq_tau_sub_7,rq_tau_sub_8,  rq_tau_w_c, covariate.labels = ROW_LABS)  ## need to change 8 for str


## Tau = 0.75
tau <- 0.75

for (i in 1:8){
  print(paste(i," ",subset_names[i]," ",sep=""))
  DATA   <- get(subset_names[i])
  
  if(QUANTILE == TRUE){
    STORE1 <- rq( EQUATION_C ,data = DATA , tau = tau)
    assign(paste("rq_tau_sub_", i,sep = ""), STORE1)
    
    rm(STORE1)}else{print("NO QUANTILE")}
  
  rm(DATA)  
}

rq_tau_e_c  <- rq(EQUATION_C , data = elec_comm_panel,  tau = tau)
rq_tau_w_c  <- rq(EQUATION_C,  data = water_comm_panel, tau = tau)

stargazer(rq_tau_sub_1,rq_tau_sub_2,rq_tau_sub_3,rq_tau_sub_4,  rq_tau_e_c, covariate.labels = ROW_LABS)
stargazer(rq_tau_sub_5,rq_tau_sub_6,rq_tau_sub_7,rq_tau_sub_8,  rq_tau_w_c, covariate.labels = ROW_LABS)


## Tau = 0.25
tau <- 0.25

for (i in 1:8){
  print(paste(i," ",subset_names[i]," ",sep=""))
  DATA   <- get(subset_names[i])
  
  if(QUANTILE == TRUE){
    STORE1 <- rq( EQUATION_C ,data = DATA , tau = tau)
    assign(paste("rq_tau_sub_", i,sep = ""), STORE1)
    
    # print(summary(STORE1, se = "iid"))
    rm(STORE1)}else{print("NO QUANTILE")}
  
  rm(DATA)  
}

rq_tau_e_c  <- rq(EQUATION_C , data = elec_comm_panel,  tau = tau )
rq_tau_w_c  <- rq(EQUATION_C,  data = water_comm_panel, tau = tau )

stargazer(rq_tau_sub_1,rq_tau_sub_2,rq_tau_sub_3,rq_tau_sub_4,  rq_tau_e_c, covariate.labels = ROW_LABS)
stargazer(rq_tau_sub_5,rq_tau_sub_6,rq_tau_sub_7,rq_tau_sub_8,  rq_tau_w_c, covariate.labels = ROW_LABS) 
